Problema 1

Dada una matriz \(4 \times 4\) de variables aleatorias Bernoulli denotada por \(\displaystyle \left[X_{ij}\right]\), sean \(N(X)\) el número total de éxitos en \(X\) y \(D(X)\) el total de vecinos (horizontales o verticales) de dichos éxitos que difieren.

Asimismo, supongamos que la distribución de \(X\), \(\pi(X)\), es proporcional a

\[ \pi(X) \propto p^{N(X)} (1 - p)^{16 - N(X)} \exp \left( -\lambda D(X) \right). \]

Si \(\lambda = 0\), las variables son independientes e idénticamente distribuidas \(\text{Bernoulli}(p)\).

Hay \(2^{16}\) posibles estados, uno por cada valor posible de \(X\). Usaremos el método Metropolis-Hastings con los siguientes kérneles de transición:

  1. \(q_1\) tal que cada transición es igualmente plausible con probabilidad \(1/2^{16}\). Es decir, el siguiente estado candidato para \(X\) es un vector de \(16\) \(\text{Bernoulli}(p)\) independientes e idénticamente distribuidas.

  2. \(q_2\) tal que se elige una de las \(16\) entradas en \(X\) con probabilidad \(1/16\), y luego se determina el valor de la celda a ser \(0\) o \(1\) con probabilidad \(0.5\) en cada caso. Entonces, a lo más puede cambiar un solo elemento de \(X\) en cada transición.

Ambas \(q_1\) y \(q_2\) son simétricas, irreducibles y tienen diagonales positivas. Es fácil ver que la primera se mueve más rápido que la segunda.

Estamos interesados en la probabilidad de que todos los elementos de la diagonal sean \(1\). Al usar \(q_1\) y \(q_2\), estimaremos dicha probabilidad para los valores \(\lambda \in \{0, 1, 3\}\) y \(p \in \{0.5, 0.8\}\).

En primer lugar, modelamos el generador de nuestra matriz de inicio y las funciones necesarias para nuestras simulaciones — \(N(X)\), \(D(X)\), \(\pi(X)\), \(q_1\) y \(q_2\).

init <- function(m, n) {
  X <- matrix(rbinom(m * n, 1, runif(1)), m, n)
  return(X)
}

N <- function(X) {
  sum(X == 1)
}

D <- function(X) {
  n <- 0
  
  error <- function(x) {
    return(1)
  }
  
  for (i in 1:dim(X)[1]) {
    for (j in 1:dim(X)[2]) {
      up     <- tryCatch(X[i - 1, j], error = error)
      down   <- tryCatch(X[i + 1, j], error = error)
      left   <- tryCatch(X[i, j - 1], error = error)
      right  <- tryCatch(X[i, j + 1], error = error)
      centre <- X[i, j]
      
      ifelse((centre == 1) & (up == 0), n <- n + 1, n)
      ifelse((centre == 1) & (down == 0), n <- n + 1, n)
      ifelse((centre == 1) & (left == 0), n <- n + 1, n)
      ifelse((centre == 1) & (right == 0), n <- n + 1, n)
    }
  }
  
  return(n)
}

dist <- function(X, p, lambda) {
  m <- dim(X)[1]
  n <- dim(X)[2]
  
  (p ^ N(X)) * ((1 - p) ^ ((m * n) - N(X))) * exp(-lambda * D(X))
}

q.1 <- function(X) {
  m <- dim(X)[1]
  n <- dim(X)[2]
  Y <- matrix(rbinom(m * n, 1, p), m, n)
  return(Y)
}

q.2 <- function(X) {
  m <- dim(X)[1]
  n <- dim(X)[2]
  a <- sample(m, 1)
  b <- sample(n, 1)
  x <- rbinom(1, 1, 0.5)
  X[a, b] <- x
  return(X)
}

Nótese que, dado que \(q_1\) y \(q_2\) son simétricas, podemos definir nuestra razón de Hastings como sigue:

\[ \alpha (x, y) = \min \left\{1, \frac{\pi(y)}{\pi(x)}\right\}. \]

simulation <- function(iterations, X, p, lambda, q) {
  m <- dim(X)[1]
  n <- dim(X)[2]
  chain <- matrix(NA, length(c(X)), iterations)
  chain[, 1] <- c(X)
  
  for (t in 2:(iterations)) {
    x <- chain[, t - 1]
    sample.x <- matrix(x, m, n)
    y <- c(do.call(q, list(sample.x)))
    sample.y <- matrix(y, m, n)
    
    U <- runif(1)
    MHR <- dist(sample.y, p, lambda) / dist(sample.x, p, lambda)
    alpha <- min(1, MHR)
    
    if (U <= alpha) {
      chain[, t] <- y
    }
    else {
      chain[, t] <- x
    }
    
    diag.count <- sum(diag(matrix(chain[, t], m, n)) == 1)
    identity <- sum(diag(diag(nrow = m, ncol = n)) == 1)
    
    if (diag.count == identity) {
      break
    }
  }
  
  keep.cols <- which(apply(!is.na(chain), 2, all))
  chain <- chain[, keep.cols]
  no.iter <- ncol(chain) # periodo de calentamiento
  final.matrix <- matrix(chain[, no.iter], m, n)
  return(list(chain, final.matrix, no.iter))
}

Únicamente usaremos la semilla para generar nuestra matriz inicial. Para estimar las probabilidades de tener \(1\)’s en la diagonal, tomaremos la media sobre \(k\) simulaciones del periodo en que se alcanzan dichos \(1\)’s en la diagonal.

k <- 100
m <- n <- 4
p <- c(0.5, 0.8)
lambda <- c(0, 1, 3)

set.seed(181121)
X <- init(m, n)

# Usando q1
for (i in p) {
  for (j in lambda) {
    prob <- mean(replicate(k, 1 - simulation(k, X, i, j, q.1)[[3]] / k))
    cat("( p =", i, ", lambda =", j, ")", "probability:", prob, "\n")
  }
}
## ( p = 0.5 , lambda = 0 ) probability: 0.9313 
## ( p = 0.5 , lambda = 1 ) probability: 0.6879 
## ( p = 0.5 , lambda = 3 ) probability: 0.5995 
## ( p = 0.8 , lambda = 0 ) probability: 0.8957 
## ( p = 0.8 , lambda = 1 ) probability: 0.788 
## ( p = 0.8 , lambda = 3 ) probability: 0.7285
# Usando q2
for (i in p) {
  for (j in lambda) {
    prob <- mean(replicate(k, 1 - simulation(k, X, i, j, q.2)[[3]] / k))
    cat("( p =", i, ", lambda =", j, ")", "probability:", prob, "\n")
  }
}
## ( p = 0.5 , lambda = 0 ) probability: 0.1599 
## ( p = 0.5 , lambda = 1 ) probability: 0.0131 
## ( p = 0.5 , lambda = 3 ) probability: 0 
## ( p = 0.8 , lambda = 0 ) probability: 0.2956 
## ( p = 0.8 , lambda = 1 ) probability: 0.1749 
## ( p = 0.8 , lambda = 3 ) probability: 0.0025

Nótese que \(q_1\) es muchísimo más eficiente que \(q_2\), pues alcanza el periodo de calentamiento mucho antes y, por lo tanto, las probabilidades de tener \(1\)’s en la diagonal respecto al número de simulaciones es mucho mayor. Asimismo, es fácil ver que una \(\lambda\) más grande alenta los procesos iterativos dada una mayor dependencia entre las variables.

Para comparar las cadenas resultantes con \(q_1\) y \(q_2\), fijemos \(p = 0.5\) y \(\lambda = 0\) para eficientar los cálculos internos.

set.seed(181121)
chain.1 <- simulation(1000, X, 0.5, 0, q.1)
set.seed(181121)
chain.2 <- simulation(1000, X, 0.5, 0, q.2)

De esta forma, podemos observar nuestras funciones de autocorrelación para nuestros elementos de interés.

Finalmente, dada nuestra matriz inicial

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    1
## [2,]    0    0    0    0
## [3,]    0    0    0    0
## [4,]    1    1    1    0

podemos ver que, con \(q_1\) y \(q_2\), nuestras matrices finales son

##      [,1] [,2] [,3] [,4]
## [1,]    1    0    1    1
## [2,]    0    1    1    1
## [3,]    0    1    1    0
## [4,]    1    1    1    1

y

##      [,1] [,2] [,3] [,4]
## [1,]    1    1    0    0
## [2,]    1    1    0    0
## [3,]    0    1    1    0
## [4,]    1    1    0    1

Mismas que se alcanzan en las iteraciones \(4\) y \(185\) (periodos de calentamiento), respectivamente. Además, es fácil observar que, la autocorrelación con \(q_1\) es menor y, por lo tanto, más estable y capaz de alcanzar la estacionariedad más rápido.

Problema 2

Queremos investigar el desempeño del algoritmo de Metropolis-Hastings usando una caminata aleatoria cuando la distribución objetivo es una mezcla de dos densidades normales bivariadas para \(\theta = (\theta_1, \theta_2)\)

\[ \pi(\theta) = 0.7 \mathscr{N} (\theta \mid \mu_1,\ \Sigma_1) + 0.3 \mathscr{N} (\theta \mid \mu_2,\ \Sigma_2) \]

donde

\[ \mu_1 = \begin{pmatrix} 4 \\ 5 \end{pmatrix},\ \mu_2 = \begin{pmatrix} 0.7 \\ 3.5 \end{pmatrix},\ \Sigma_1 = \begin{pmatrix} 1 & 0.7 \\ 0.7 & 1 \end{pmatrix},\ \Sigma_2 = \begin{pmatrix} 1 & -0.7 \\ -0.7 & 1 \end{pmatrix}. \]

Supongamos que no podemos muestrear \(\pi\) directamente y, en su lugar, usamos Metropolis-Hastings con

\[ q(y \mid x) = \mathscr{N}(y \mid x, v I_2) \]

donde \(v\) es un parámetro de ajuste.

mu.1 <- c(4, 5)
mu.2 <- c(0.7, 3.5)
sigma.1 <- matrix(c(1, 0.7, 0.7, 1), 2, 2)
sigma.2 <- matrix(c(1, -0.7, -0.7, 1), 2, 2)
I.2 <- diag(2)

mix.norm <- function(x) {
  0.7 * dmvnorm(x, mean = mu.1, sigma = sigma.1) +
    0.3 * dmvnorm(x, mean = mu.2, sigma = sigma.2)
}

q <- function(a, b, c) {
  dmvnorm(a, mean = b, sigma = (c * I.2))
}

# Simulación
simulation <- function(n, init, v) {
  X <- matrix(c(NA, NA), 2, n)
  X[, 1] <- init
  count <- 0
  
  for (t in 1:(n - 1)) {
    x <- X[, t]
    y <- rmvnorm(1, mean = x, sigma = (v * I.2))
    U <- runif(1)
    MHR <- (mix.norm(y) * q(x, y, v)) /
      (mix.norm(x) * q(y, x, v))
    alpha <- min(1, MHR)
    
    if (U <= alpha) {
      X[, t + 1] <- y
      count <- count + 1
    }
    else {
      X[, t + 1] <- x
      count <- count
    }
  }
  
  prob <- count / n
  X <- data.frame(x = X[1, ], y = X[2, ], n = 1:n, prob = prob, v = v)
  
  return(X)
}

# Visualización
visualize <- function(X) {
  colours <-
    c("Punto inicial" = "#1e40ca",
      "Media posterior" = "#00a2ed")
  
  p.1 <- ggplot(X) +
    geom_path(aes(x = x, y = y), size = 0.1, alpha = 0.3) +
    geom_point(aes(x = x, y = y), size = 0.1) +
    geom_point(aes(x = x[1],
                   y = y[1],
                   colour = "Punto inicial"), size = 2) +
    geom_point(aes(
      x = mean(x),
      y = mean(y),
      colour = "Media posterior"
    ), size = 2) +
    labs(
      title = NULL,
      x = expression(theta[1]),
      y = expression(theta[2]),
      caption = paste("Tasa de aceptación:", X$prob[1],
                      "\nv =", X$v[1])
    ) +
    scale_colour_manual(name = NULL, values = colours) +
    theme(legend.position = "top")
  
  p.2 <- ggplot(X) +
    geom_line(aes(x = n, y = x),
              size = 0.1) +
    labs(title = NULL,
         x = "Iteración",
         y = expression(theta[1]))
  
  p.3 <- ggplot(X) +
    geom_line(aes(x = n, y = y),
              size = 0.1) +
    labs(title = NULL,
         x = "Iteración",
         y = expression(theta[2]))
  
  p.4 <- ggplot(X, aes(x = x, y = ..density..)) +
    geom_histogram(binwidth = 0.5,
                   fill = "#1e40ca",
                   alpha = 0.6) +
    geom_density() +
    labs(title = NULL,
         x = expression(theta[1]))
  
  p.5 <- ggplot(X, aes(x = y, y = ..density..)) +
    geom_histogram(binwidth = 0.5,
                   fill = "#1e40ca",
                   alpha = 0.6) +
    geom_density() +
    labs(title = NULL,
         x = expression(theta[2]))
  
  plots <- list(p.1, p.2, p.3, p.4, p.5)
  return(plots)
}

A

Al momento de otorgar valores pequeños o grandes a \(v\) en \(q(y \mid x)\), el efecto de dicha \(v\) está directamente relacionado con nuestra razón de Hastings, pues esta se escala a varianza \(v\), reduciendo o aumentando la tasa de aceptación. Una \(v > 1\) implica, dada la definición de \(\alpha(x, y)\), una tasa de aceptación menor. Asimismo, una \(v < 1\) implica, naturalmente, una tasa de aceptación mayor.

Nótese lo anteriormente mencionado en las siguientes caminatas aleatorias para \(v = 0.01\) y \(v = 100\):

set.seed(181121)
visualize(simulation(1000, c(0, 0), 0.01))[[1]]

set.seed(181121)
visualize(simulation(1000, c(0, 0), 100))[[1]]

B

Función de autocorrelación de \(\theta_2\):

set.seed(181121)
Y <- simulation(1000, c(0, 0), 1)
acf(Y$y, main = expression(theta[2]))

C

Para \(5,000\) extracciones, notemos que nuestra tasa de aceptación es de \(0.4788\).

Gráficamente:

Problema 3

Sea \(\displaystyle \theta = \int_0^{\pi/3} \sin(t)dt\), usaremos Monte Carlo para calcular un estimador \(\hat{\theta}\).

Monte Carlo crudo

a <- 0
b <- pi / 3
n <- 100000

set.seed(181121)
sim.u     <- runif(n, a, b)
sim.raw   <- sin(sim.u)
theta.raw <- cumsum(sim.raw) / (1:n)
theta.raw <- (b - a) * theta.raw
theta.1   <- tail(theta.raw, 1)
var.1     <- var(sim.raw)

Nuestro estimado es: \(\hat{\theta}_1 = 0.5003\).

Variables antitéticas

sim.v     <- a + (b - sim.u)
sim.ant   <- (sin(sim.u) + sin(sim.v)) / 2
theta.ant <- cumsum(sim.ant) / (1:n)
theta.ant <- (b - a) * theta.ant
theta.2   <- tail(theta.ant, 1)
var.2     <- var(sim.ant)

Nuestro estimado es: \(\hat{\theta}_2 = 0.4999\). Asimismo, con variables antitéticas, la varianza se reduce en un \(99.3869\%\).

Variables de control

Sabemos \(\sin(t) = \cos(\pi / 2 - t)\). Definamos \(x = \pi / 2 - t\) como nuestra variable de control.

k <- 1000

# Piloto
set.seed(181121)
u <- runif(k, a, b)
x <- pi / 2 - u
y <- cos(x)
alpha <- -lm(y ~ x)$coeff[2]

# Simulación
set.seed(181121)
u <- runif(n, a, b)
x <- pi / 2 - u
y <- cos(x)

sim.cv   <- y + alpha * (x - mean(x))
theta.cv <- cumsum(sim.cv) / (1:n)
theta.cv <- (b - a) * theta.cv
theta.3  <- tail(theta.cv, 1)
var.3    <- var(sim.cv)

Nuestro estimado es: \(\hat{\theta}_3 = 0.5003\). De tal forma que, con nuestra variable de control, la varianza se reduce en un \(99.3659\%\). Finalmente, nuestro estimador \(\hat{\theta}_3\) reduce su varianza en un \(3.4153\%\) respecto a \(\hat{\theta}_2\).

Gráficamente,

Problema 4

Dados los siguientes números normales con \(\mu = 3\) y \(\sigma = 1\)

obs <- c(4.59, 4.153, 2.46, 2.732, 2.973)

la función de distribución empírica \(F_5(x)\) está dada por:

\[ F_5(x) = \begin{cases} 0, & \text{si}\ \ x_{(1)} > x \\ 1/5, & \text{si}\ \ x_{(1)} \leq x < x_{(2)} \\ 2/5, & \text{si}\ \ x_{(2)} \leq x < x_{(3)} \\ 3/5, & \text{si}\ \ x_{(3)} \leq x < x_{(4)} \\ 4/5, & \text{si}\ \ x_{(4)} \leq x < x_{(5)} \\ 1, & \text{si}\ \ x_{(5)} \leq x \end{cases} \]

Asimismo,

\[ \begin{align} P(F_5(3) \leq 2/5) &= P(F_5(3) = 0) + P(F_5(3) = 1/5) + P(F_5(3) = 2/5) \\ &= {5 \choose 0} F(3)^0 \left( 1 - F(3) \right)^{5 - 0} \\ &+ {5 \choose 1} F(3)^1 \left( 1 - F(3) \right)^{5 - 1} \\ &+ {5 \choose 2} F(3)^2 \left( 1 - F(3) \right)^{5 - 2} \end{align} \]

donde \(F(x)\) es una normal con \(\mu = 3\) y \(\sigma = 1\).

Así pues,

mean <- 3
sd <- 1

prob <- 0

for (i in 0:2) {
  p <- choose(5, i) *
    (pnorm(3, mean, sd)) ^ (i) *
    (1 - pnorm(3, mean, sd)) ^ (5 - i)
  
  prob <- prob + p
}

prob
## [1] 0.5

De forma analítica, es fácil ver que \(F(3) = \Phi(0) = 0.5\) donde \(\Phi(\cdot)\) es una normal estándar. De igual modo, nótese que \(F_5(x) \leq 2/5\) si y solo si \(x < x_{(3)}\) donde \(x_{(3)}\) es la mediana de nuestra muestra. Por lo tanto, podemos concluir que \(P(F_5(3) \leq 2/5) = 0.5\).